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ABSTRACT 

It is well known that the classical gravitational two body problem can be transformed 
into a spherical harmonic oscillator by reguiarization. We find that a modification 
of the reguiarization transformation has a similar result to leading order in general 
relativity. In the resulting harmonic oscillator, the leading-order relativistic perturba- 
tion is formally a negative centrifugal force. The net centrifugal force changes sign 
at three Schwarzschild radii, which interestingly mimics the innermost stable cir- 
cular orbit (ISCO) of the full Schwarzschild problem. Transforming the harmonic- 
oscillator solution back to spatial coordinates yields, for both timelike and null weak- 
field Schwarzschild geodesies, a solution for t, r, (f> in terms of elementary functions 
of a variable that can be interpreted as a generalized eccentric anomaly. The text- 
book expressions for relativistic precession and light deflection are easily recovered. 
We suggest how this solution could be combined with additional perturbations into 
numerical methods suitable for applications such as relativistic accretion or dynamics 
of the Galactic-centre stars. 
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1 INTRODUCTION 

When Schwarzschild geodesies appear in classic tests of general relativity, the important result is an integral over the geo desic: 
orbit al precession or deflection of light. Similary, in modern tests of relativity involving binary pulsars (for a review, see WilJ 
2006l l the observable effects are also cumulative over many orbits. 

In the case of the recently-discovered S-stars near the Galactic Centre, the cumulative effects of relativity are no lon ger the 
principal quantity of interest. The highly eccentric examples from the S-stars ( Ghez et al. 20081 : Gillessen et al. 20091 ) . which 
experience a range of gravitational regimes, motivate an interest in tracking relativistic effects as they vary al ong an orbit 



In particular, some recent work has drawn attenti on to relativistic effects on redshifts near pericentre passage ijZucker et al. 



20061 : iKannan fc s"ahall2009l : lAngelil fc Sahalboiol ). These effects can be calculated numerically and some of them also by 
post-Newtonian perturbation theory, but a simpler method is desirable. 

Such a method is suggested by Levi-Civita or Kustaanheimo-Stiefel (LC or KS) reguiarization, which are transformations 
of the classical gravitational two-bo dy problem t o an equivalent harmonic oscillator. This type of reguiarization was o rigi- 
nally introduced in two dime nsions ( Levi-Civitalll920.1 and much later extended to three dimensions (Kustaanheima ,19641 : 
Kustaanheimo fc Stiefellll965). KS reguiarization has an extensive literatur e, including applications to A'^-body simulations 
I Aarseth fc Zarelll974al lbl: ljernigan fc Porteil[l989l : iMikkola fc Aarseth|[l993 ) . The classical result suggests that the LC or KS 
reguiarization could be used to transform the general relativistic problem into a perturbed harmonic-oscillator. We find even 
better: a modification of the LC /KS transforma tion acting on the geodesies of the leading order Schwarzschild metric in the 
isotropic or harmonic gauge (cf. IWeinberg||l972l . Section 8.2), 
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yields an unperturbed circular/spherical harmonic oscillator. As a result, the solution is analytic. The difference from the 
classical case is a negative centrifugal- force term in the transformed space. This terms encodes the leading-order effects of 
precession, deflection of light, and the innermost stable orbit. 

Because orbits in the Schwarzschild spacetime do not leave the orbital plane, in this paper we mainly consider the 
two-dimensional or LC case. The three-dimensional or KS case is similar, but algebraically more complicated, as it involves 
introducing a fourth spatial dimension. 



2 TRANSFORMATION OF THE GEODESIC EQUATIONS 

Since regularization is formulated in the language of Hamiltonians, we begin by exp ressing the geodesic equations in Hamil- 
tonian form. A convenient expression for the Ifamiltonian is (e.g.. Equation 25.10 of lMisner et~al]|l973h 

H = ^S^VmP- (2) 
Considering the metric fl]), we have 



1-^)^, (3) 



'' = -2[' + — + —)^^ +V-—J 2 

This retains the leading-order Newtonian corrections, which are 0{M/r) in the spatial part and 0{M'^ /r^) in the temporal 
part, and neglects higher orders. Note that the Ifamiltonian ((3]) has four degrees of freedom: time t being a coordinate and pt 
its conjugate momentum, while the affine parameter A is the independent variable. We assume units with G/(? — 1. 
The Hamiltonian being independent of t, it follows that pt is a constant, and Hamilton's equation for t reads 

dH dt / 2M 2M'^\ 

where A denotes the afhne parameter. On choosing pt = —1 (which we are free to do, as this amounts to choosing units for 
A) and discarding a constant we arrive at the Ifamiltonian 

1 / 2M 2A/2\ / 2M\p^ 
'' = -2[—^—)^[^-—)-2- 

As is well known from the spherical symmetry of the Schwarzschild spacetime, geodesies are confined to a plane. Without 
loss of generality we can choose a two dimensional coordinate system in the x, y plane. (Below, at the end of this section, we 
briefly indicate the procedure without this simplification.) 

We now apply a regularization transformation, fn the x, y plane, we introduce two new coordinates, given by the real 
and imaginary parts of the complex number 



Q = \/x + iy- (6) 

Hence r — \Q\'^ ■ The conjugate momentum components are given by the real and imaginary parts of a complex number P, 
which satisfies 

P.+ip, = ^. (7) 
On multiplying the above equation by its complex conjugate we arrive at the transformation 

p|2 

The P, Q are know n as Levi-Civita (or LC) variables. A proof that they are indeed canonical appears in several sources (for 
example lSahalboogl ) and we do not repeat it here. Transforming ((5]) and rearranging yield the Hamiltonian in the LC variables 

1 /2M 2M'\ f 1 2Af\ \P\^ 

To complete the regularization we invoke a Poincare time transformation which involves introducing a scaled time variable 
related to the affine parameter by 

d\^g{P,Q)ds, (10) 
where g{P, Q) can be any function, and define a new Hamiltonian 

r = g{H-E). (11) 
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The r Hamiltonian preserves the Hamihonian form of the equations of motion, provided the constant E is the initial 
value of the original Hamiltonian. We choose 

,= ('^_MV' = (12) 

which approximates to \Q\^ + 2M for large \Q\. The time-transformed Hamiltonian takes the form 

lp|2 a.r2 

r = H_L _ i5jQ|2 - - M(i + 2i5) + o (lor^) . (13) 

We remark that the classical case uses g — \Q\'^, and only the first two terms in (|13p are present in the F Hamiltonian, 
up to a constant. However, as is evident from (|13p . the weak-field Schwarzschild case also yields a harmonic oscillator under 
a regularization transformation. The essential difference is in the extra 1/|Q|^ term which encodes the leading order effects 
of general relativity by altering the classical angular momentum in (I13p . Thus, the key modification from the classical case 
which allows the formulation of an analytically solvable, relativistic Hamiltonian is the choice of g{P, Q) in equation (|12|l . 

If we do not restrict the coordinates to the plane, the LC transformation must be replaced with a Kustaanheimo-Steifel 
transformation. Although the transformation itself is far more complicated (involving four spatial dimensions) the time 
transformation and F Hamiltonian (|13|) remain the same, except that |P| and \Q\ are lengths in four Euclidean dimension s. 
This is easily seen on comparing the above with KS regularization of the Kepler problem (see, e.g.. Section 5 of Saha 20091 ). 



3 SOLUTIONS 

Continuing in two dimensions, it is possible to put the F Hamiltonian in a more recognizable form by transforming to polar 
coordinates {Qr,Q4,, Pr, P,/,), 

p2 (Pj-6(2Af)2) „ , 

r = ^ + ^ ' - EQl - M{1 + 2E) + {Q~^) , (14) 

which for negative values of E is the Hamiltonian for a circular classical harmonic oscillator with squared angular momentum 
decreased by 6(2M)^ from the equivalent Kepler problem. 

We remark that the F Hamiltonian appears to depend only on terms up to 1/r, and it would seem that by leaving out 
terms of order we have omitted relativistic effects. However, this is not the case, since the time equation now hides 

a factor of r. Thus, one should read terms of order IQI " as terms of order r"'"''^^^' for n = 0, 1, 2 . . . This provides some 
insight as to why the Hamiltonian now has soluble equations of motion: we have pushed a factor of r into the time equation. 

We introduce the constants 

M' = M{l + 2E) and P'l = P^ - 6{2Mf . (15) 
Dropping higher order terms, the Hamiltonian becomes, 

r = f + g-i^o?-M', (16) 

which is identical to the transformed Hamiltonian of a particle with angular momentum P^ in a central force potential of the 
form V{r) — M' /r, or equivalently to the Kepler problem where A/ is replaced by Af' and P| is replaced by P'^. 
Hamilton's equations of motion are 

"ds~ ^ T 



,/2 

2EQr + -^-^ (17) 



dPr _ , p'; 1 

— - 2EQ. + —^ 



dQcf, _ P^ 1 

while P^ is a constant. On combining with the time transformation (|10|l . the equation Q for the time coordinate becomes 



, . , ^ 2Ar^ \ ds 
dt = [Ql + 2M+—^]———, 



Q2 / l-2Af/Q2 



AM + Ql + G'-^+Oi 



.A/2 



where we have again used the approximation ^ 2Af 



ds, (18) 
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We remark that = "ipif,, which is the weak- field relationship between LC angular momentum and angular momentum 
in standard coordinates. This can be seen from a comparison of the third of equations (|17p and the corresponding Hamilton 
equation applied to ((5| in the time variable s. 



3.1 Bound orbits 

For negative values of E the solutions are bound orbits. We define a 'classical' and a 'relativistic' semi-major axis 
M , M' 

the eccentric anomaly, 
and the eccentricity 



1 + 



P'lE 



(19) 



(20) 



(21) 



2Af'2 

in the standard way. Solving equation (|18|) gives an implicit equation for the coordinate time t in terms of /3 and solving 
equations (|17|l via quadrature gives the following solutions for the LC variables in terms of elementary functions of jS. 

QliP) = a'[l-ecos/?] 



-tan 



1 + e 
tan 

1 — e 



(22) 



t = 



M' 



AM 



■■13 + 



6(2A/)2 



Tmin- The LC radial momentum Pr(/3) is then generated from the 



where we have chosen Q4,{0) — 0, and Qr{0) ~ a;'(l — 
first of equations (|17|l and, with the above choice of Q^{0), the initial radial momentum vanishes. 

Note that in terms of the classical semi-major axis of the Kepler problem the quantity a — a 
show that equations (|22|) reduce to the Kepler LC equations of motion when M is small. 



M. So it is simple to 



3.2 Unbound trajectories 

For the unbound case the Hamiltonian and thus the equations of motion only change by the sign of E which appears in a', 
e, and /3. With this substitution the unbound solutions become 

QliP) = a' [e cosh /3 - 1] 



Qm = ^tan-^ 



e+ 1 



itanh(0 



(23) 



t = 



— [esmh/3 - /3] + -^=13 + ^ ^ ' Q^, 



where we have chosen Ql{0) = a' (e — 1), which in the unbound case is the point of closest approach. 



3.3 Light rays 

Null geodesies are the solutions to Hamilton's equations when H in equation ([3]) is set to zero. Since we discarded a constant 
of —1/2 in the derivation of the LC Hamiltonian, the null solutions can be found by assigning E in the unbound case the 
value of 1/2. This amounts to redefining the constants a', /3, and e in equations (|23p such that the null equations of motion 
become 



QI{s) 
Q4,{s) 

t 



2M [e„ cosh s - 1] 



■ tan 



e„ -f 1 



tanh 



2M [e„ sinh s - s] + 4M s 



&(2Mf 



(24) 
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with 
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and Q^(o) = 2M(e„ - 1). 



4 PROPERTIES OF THE SOLUTION 

Now that we have derived the bound, unbound, and nuU equtions of motion we show that, to first order, they reproduce 
the predictions for geodesies in a Schwarzschild spacetime, namely those of orbital precession, the deflection of light, and the 
innermost stable circular orbit (ISCO). 



4.1 Orbital precession 



The prefactor of the middle line in (|22p automatically gives the precession rate of orbits and tells us that the precession is due 
to the —6{2M)^ perturbation of the classical squared angular momentum. This is equivalent to the conventional interpretation 
of precession being caused by an additional centrifugal force term. Substituting for PL we have 



PL 



1 



,{2Mf 



1 -1/2 



Pi 



Since P| oc Qf, then P| >• (2M)^ and we may write 



1-6 



(2M)2 



1 -1/2 



PI 



{2Mf 

4> 



(26) 



(27) 



We note that, due to the complex square root nature of the LC transformation, — ^(j). Thus as the solution evolves through 
one period will increase by 



(2Mf 
TV + Stt — -^rrj — rad 



(28) 



giving an orbital precession rate of 
{2Mf rad 



AQd, — Stt- 



p2 



orbit 



(29) 



Converting back to non-LC coordinates and expressing the previous equation in terms of the semi-latus rectum of the orbit 
Q = a (1 — e), where a is the classical semi-major axis, the precession rate becomes 

MG rad. 



A(f: = 6tt 



c^a orbit' 



(30) 



as in the conventional treatment (e.g., Weinberg 19721 ). 



4.2 Deflection of light 



We may similarly derive the deflection angle of a light ray. Imagine that the light ray starts infinitely far from the Schwarzschild 
mass at i = — oo, approaches the point of closest approach at t = 0, and continues on to t = oo. If there is no defiection 
will sweep out n/2 radians since — ^0. If there is defiection the total difference in Q^{~oo) and Q^{oo) will be greater 
than ^. To determine the deflection angle we compute. 



AQ<^ = 2[Q^(oo)-Q^(0)] 



TT _i 

— = 2-^ tan 



+ 1 







(31) 



where is the null solution for Q^. Substituting for P'^ and the expression for e„ in equation (|25|l . invoking the binomial 
approximation for P'^ 3> M'^ and Pj ^ A/^, and neglecting terms of order and greater we find 



AQ^, = 2 (l + O (M^)) tan"^ 



1 



— + O (M^) 



+ 



Using the expansion 

tan"^ {1 + x) = j + ]^x + (x^) 

we obtain 



(32) 



(33) 
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+ 2 



2M 



0{M'^ 



TV 

2' 



(34) 



Converting back to non-LC coordinates we gain a factor of 1/2 since the LC angular momentum is twice the non-LC 
angular momentum p'^. We may treat p'^ as the corrected (relativistic) angular momentum which is equal to the impact 
parameter b for photons (with the speed of light set to unity). Thus we arrive at the expression, to first order in M, for the 
angle by which light is deflected due to a spherically symmetric mass distribution, 

4MG AUG 



Where to this order the impact parameter can be replaced by ro, the point of closest approach. 



(35) 



4.3 Innermost stable orbits 

Due to the weak-field approximation used to derive the geodesic equations of motion, the solutions do not exhibit an event 
horizon. Remarkably though, the solutions do reproduce the phenomenon of an innermost stable circular orbit. This is due 
to the centrifugal term Q{2M)'^ /Q^ in the LC Hamiltonian p4p . 

From the first two of Hamilton's equations (|17[) for a bound orbit we find 



d^Qr \E\ 



ds^ 2 



1 - 



p/2 

'>\E\Q* 



(36) 



and notice that equation (|36|) becomes a one dimensional harmonic oscillator in the radial coordinate (i.e., the condition for 
radial free fall, or equivalently e — >■ 1) when P'^ = 0. Writing out P^, the condition becomes 

P| = 6(2M)^ (37) 

Dividing by \E\ and using the definitions of the classical semi-major (a) and semi-minor (6) axes of the orbit this condition 
becomes 

= 6aM. (38) 

Which in the case of a circular orbit (a = 6) is precisely the ISCO predicted by the conventional analysis of timelike geodesies 
in a Schwarzschild metric. Note that although this relation has been derived for the classical semi-major axis one can convert to 
the standard radial coordiante of the Schwarzschild metric and find an identical result as follows. First convert to the relativistic 
semi-major axis, a' — a — M which, for the above ISCO, gives a' = 5M. Then recall that the standard Schwarzschild radial 
coordinate R is given in terms of the isotropic radial coordinate by (see e.g. Weinberg 1972 section 8.2) 

where the approximation is made to be consistent with the derivation of the equations of motion. For a' = r = 5M this gives 
R — 6M as expected. 

To interpret the ISCO derived here we observe from equation (|36|) that the condition for zero radial acceleration is 



PJ-6(2M)2 

Now consider a classical orbit with a fixed P| and imagine 'turning on' relativity adiabatically keeping P| constant. As 
relativity is turned on the negative centrifugal force term increases and effectively reduces P'^. If P^ > 6(2M)^ to begin with, 
then the classical orbit will shrink by the appropriate amount in the presence of relativistic effects. But if we originally had 
^ 6(2M)'^, then as relativity is turned on equation H40|) shows that would shrink to implying that all orbits for which 
P| < 6{2Mf are unstable. 

Also note that the standard coordinate singularity at P = 2M is mapped to r = AI in the aproximate isoptropic 
coordinates. This yields the following interpretation for the the LC equations of motion. Substituting a' = a ~ M in the first 
of equations (|22[) for the LC radial coordinate and restricting to the circular case where e = we find 

Qlil3) + M = a, (41) 

which has the interpretation that the LC radial coordinate is measured not from the origin of the corresponding Kepler 
problem but from the event horizon predicted by the Schwarzschild metric. This does not carry the interpretation of a horizon 
since when e 7^ orbits can still come arbitrarily close to r = and continue outside of r = M. However, it does add an 
interesting interpretation to the regularization transformation. For circular orbits the regularization essentially cuts out the 
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area inside the event horizon and stitches it back together mapping a circle to the origin. For elhptical orbits one can imagine 
an analogous interpretation. 

Since these equations of motion predict the correct ISCO and exhibit special behavior at the Schwarzschild event horizon, 
perhaps they could be particularly useful for simple approximate modelling of relativistic accretion discs. 



4.4 Error terms 



Working backwards from the approximate Ifamiltonian in the time variable s (equation I13p by solving equation (|lip for H 
we can recover the error terms in our approximate Hamiltonian. Solving 

\Q\' 



r = 



:{H-E) 



l-2M/|Qi2 

with no further approximations and converting to non-LC coordinates yields what we may call a surrogate Hamiltonian 



(42) 



HsUTT 



1 + 



2M 2M^ f 12M^ 



V 



Pt 



I 2M \ p2 ^ 



(43) 



for which the analytic solution is exact. For timelike geodesies E is small and for null geodesies is a second order 

correction. Thus to obtain a first order expression for error in the metric components we may neglect the term containing E. 

If it is bothersome that the Hamiltonian in equation (|43p is dependent on E, one can eliminate it by a further modification. 
Let us add a E/\Q\'^ term to the F Hamiltonian, 



, 4M^ 

which makes only a higher-order change to the solutions. Introducing -ffsmr by 
F' = \Q\\l + 2M/\Q\^ + 4A/VlQh (i/su„ - E) 
one finds 



H' — 



1 + 2M/r + 4M2 /r2 



2 1 
Pt + 



P_ 

2 



(44) 



(45) 



(46) 



1 + 2M/r + 4M2/r2 

which, when expanded in 1/r to the appropriate order, is identical to (|43|l less the term dependent on E. The Hamiltonian 
(|46|l is then the Hamiltonian for which equations (I22|l . (|23p . and (1241) are exact solutions. 



5 DISCUSSION 



We have derived timelike and null geodesies in the leading-order Schwarzschild metric in terms of elementary functions. The 
expressions (|22p for bound orbits and (|23p for unbound orbits, together with (|24p for light rays, are all simple generalizations 
of well-known expressions in classical celestial mechanics. The usual formulas for relativistic orbital precession and light 
defiection are easily recovered. A feature resembling the innermost stable circular orbit in the full Schwarzschild metric is also 
present. 

The technique we have used is a modification of the Levi-Civita or Kustaanheimo-Stiefel regularization transformation 
and transforms the geodesic equation into a spherical harmonic oscillator. The simplicity of the result, notwithstanding the 
non-trivial route used to derive it, hints at some underlying symmetry in the Schwarzschild problem. We speculate that it is 
somehow related to t he separability of the Hamilton- Jacobi and other equations in the Schwarzschild and Kerr metrics (cf. 
Chandrasekhai 19831 ) but have not attempted to investigate this. 

As mentioned in the Introduction, the original motivation for this work was to find useful formulas applicable to the 
highly-eccentric Galactic-centre stars, whose orbits pass through a large range of gravitational regimes. Future observations 
of these stars aiming to detect relativistic effects will require computation of relativistic effects on both stellar orbits and light 
rays at many points along an orbit, for many orbits, in order to fit the orbital parameters. The solutions in this paper allow a 
simpler, more efficient method for carrying out those computations. The analytic solutions will not be suf ficient on their own 
becau se the Galactic-centre stars also experience additional Newtonian perturbations due to local matter (IMikkola fc Merritt 
2008h , but they can be incorporated into numerical methods, specifically, generalized leapfrog i ntegrators. Such algorithms 
evol ve alternately under tw o Hamiltonians, which are integrable separately. The idea g oes back to Wi sdom fc HolmanI (1991) 
and iKinoshita et al.l (|l99lh . Some recent develop ments on adaptive step sizes appear in lEmel'vanenkol (|2007h and are applied 
to the specific problem of Galactic-centre stars in lPreto fc Sahal ( 20091 ) . We note, however, that the present work is limited 
to test particles, and hence will not be applicable for binary orbits or self-gravitating disc simulations unless a generalization 
is found. 
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Another potential application may be the use of t he solutions in relativistic disc simulations as an alternative to th e 
widely used pseudo-Newtonian potentials (see especially Paczviiski fc Wiita 1980l : Artemova et al. 19961 : Abramowicz 20091 ) . 
an advantage being that the solutions in this paper are well-defined approximations and include a more complete repertoire 
of general-relativistic effects for the same computational budget. 
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